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Abstract 

The  issue  of  high  speed  impact  and  hypervelocity  gouging  in  particular  is  of  great 
interest  to  the  Air  Force  Office  of  Scientific  Research  and  the  Air  Force  Research 
Laboratory’s  Holloman  High  Speed  Test  Track.  Rocket  sled  tests  at  the  facility 
frequently  approach  velocities  where  concern  for  gouging  development  becomes  the 
limiting  factor  to  achieving  higher  operating  velocities.  Direct  observation  of  the 
gouging  process  is  not  possible  so  computational  modeling  is  necessary  to  study  the 
phenomenon.  Since  gouging  development  is  dependent  on  the  impact  surface  conditions, 
the  method  used  to  model  material  interfaces  directly  affects  the  accuracy  of  the  solution. 

This  research  is  primarily  interested  in  the  methods  for  modeling  material 
interfaces  available  in  CTH,  a  hydrocode  developed  by  Sandia  National  Laboratories. 
Three  methods  are  available  in  CTH  to  handle  material  interfaces:  1)  materials  are  joined 
at  the  interface,  2)  a  frictionless  slide  line  is  inserted,  and  3)  a  boundary  layer  interface  is 
established.  An  axisymmetric  impact  scenario  is  used  to  explore  these  methods  and  their 
influence  on  high  energy  impact  solutions.  The  three  methods  are  also  compared  in  an 
axisymmetric  sliding  scenario. 

The  method  of  joining  materials  at  the  contact  surface  appears  well  suited  to  the 
simulation  of  high  energy  impact  events  and  hypervelocity  gouging.  In  general,  the 
frictionless  slide  line  method  introduces  numerical  instability.  Finally,  the  boundary 
layer  interface  method  is  too  limited  for  two  dimensional  applications. 
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ANALYSIS  OF  COMPUTATIONAL  METHODS  FOR  THE  TREATMENT  OF 

MATERIAL  INTERFACES 


I.  Introduction 


Background 

The  Holloman  High  Speed  Test  Track  (HHSTT)  is  a  ground-based  test  facility 
operated  by  the  US  Air  Force  for  the  testing  of  various  hardware  systems  in  simulated 
free  flight  conditions.  The  facility  conducts  tests  using  rocket  propelled  sleds  that  ride 
along  continuous  steel  rails  which  span  a  distance  of  50,988  feet.  The  rocket  sleds  are 
guided  along  the  rail  by  steel  slippers  which  maintain  a  small  gap  with  the  rail  to  allow 
for  limited  free  flight.  Upgrades  to  the  facility  have  allowed  tests  to  be  conducted  at 
extremely  high  velocities  approaching  1-2  km/s  with  plans  presently  developed  that 
would  allow  for  a  velocity  of  3 km/s.  At  these  velocities,  the  free  flight  conditions  allow 
for  the  possibility  that  the  slippers  will  impact  the  rail  due  to  roll,  pitch,  or  yaw  motions 
experienced  by  the  test  sled.  The  severity  of  the  impact  varies  but  in  certain  cases  can 
result  in  the  loss  of  the  test  sled.  Inspections  of  the  rails  after  test  events  have  revealed 
damage  in  the  form  of  gouges.  Gouges  associated  with  high  energy  impacts  are 
characterized  by  a  teardrop  shaped  shallow  depression  in  the  material  along  with  a  raised 
hump  at  the  blunt  end  of  the  gouge.  Figure  1  shows  a  diagram  of  a  typical  gouge. 
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Top  View 


Side  View 


Direction  of  Travel 


Figure  1.1  Illustration  of  a  Gouge 


A  hypothesis  for  the  hypervelocity  gouging  process  was  offered  by  Barber  and 
Bauer  [2]  which  attributed  gouge  fonnation  to  microscopic  surface  asperity  impact 
between  a  moving  slider  and  a  stationary  rail.  The  impact  is  characterized  as  a  discrete, 
localized,  and  violent  event.  The  energy  released  in  such  an  impact  would  resemble  a 
microscopic  explosion  and  produce  a  small  crater  at  the  surfaces  of  the  contacting  solids. 
The  relative  motion  of  the  slider  imparts  a  velocity  to  the  center  of  mass  of  the  explosion 
of  approximately  one-half  of  the  slider  velocity.  The  motion  of  the  center  of  mass  of  the 
explosion  as  it  expands  results  in  the  characteristic  teardrop  shaped  gouge  which  is 
formed  in  both  the  slider  and  the  rail. 


Computational  Hypervelocity  Gouging  Research 

The  hypervelocity  gouging  process,  as  it  pertains  to  the  HHSTT,  was  successfully 
simulated  by  Laird  and  Palazotto  [10]  using  the  hydrocode  CTH.  They  were  able  to 
observe  the  development  of  a  gouge  and  identify  the  mechanisms  that  lead  to  gouge 
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formation.  Two  conditions  were  reported  as  necessary  for  gouge  initiation.  First,  the 
slipper  must  penetrate  into  the  rail,  through  a  vertical  impact  or  a  horizontal  impact  with 
an  asperity.  Secondly,  the  slipper  must  maintain  sufficient  horizontal  velocity  after 
penetration  to  produce  stresses  exceeding  the  flow  stress  of  both  the  slipper  and  rail 
materials.  The  resulting  plastic  flow  takes  the  form  of  material  jets  which  continually 
grow  and  penetrate  further  into  the  two  surfaces,  initiating  a  gouge.  The  gouge  grows  in 
the  direction  of  the  slipper  velocity  but  at  one  half  of  the  magnitude  of  the  slipper 
velocity.  The  length  and  depth  of  the  resulting  gouge  are  proportional  to  the  length  of  the 
slipper  and  the  amount  of  time  required  for  the  slipper  to  move  past  the  gouge.  A  typical 
gouge  being  developed  is  shown  in  Figure  1.2. 


Figure  1 .2  Gouge  Simulation  by  Laird  [9] 


Similar  research  was  also  conducted  by  Szmerekovsky  [14]  in  which  the 
capabilities  of  three  computational  tools  to  simulate  hypervelocity  gouging  were 
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evaluated.  The  hydrocode  CTH  was  the  only  software  package  which  employed  an 
Eulerian  Finite  Volume  solution  method.  ABAQUS  and  LS-DYNA3D  were  the  other 
two  software  packages  and  both  employed  a  Lagrangian  Finite  Element  solution  method. 
The  following  1 1  criteria  were  used  by  Szmerekovsky  to  evaluate  the  three  computational 
tools  and  detennine  which  was  most  suitable  for  numerical  simulation  of  the 
hypervelocity  gouging  phenomenon  [14]. 

1 .  Mesh  capabilities  and  their  ability  to  capture  the  large  deformations 
and  high  strain  rates  of  the  gouging  phenomenon  without  excessive 
distortion. 

2.  Availability  of  equation  of  state  and  constitutive  model  capabilities. 

3.  Representation  of  failure  and  availability  of  failure  models. 

4.  Ability  of  the  solution  methods  and  material  models  to  accurately 
model  high  energy  impact  phenomena  such  as  normal  shock  waves. 

5.  Availability  of  models  for  heat  flux  and  heat  flow  due  to  deformation. 

6.  Availability  of  contact  and  material  mixing  algorithms. 

7.  Capability  to  model  conditions  (e.g.,  plastic  zone,  shear  bands,  and 
high  pressure  core)  that  lead  to  gouging. 

8.  Capability  to  model  interacting  material  “jets”  and  ejecta  associated 
with  gouging. 

9.  Availability  of  the  code  on  readily  accessible  high  performance 
computational  platfonns. 

10.  Availability  of  technical  support. 

11.  Previous  work  done  on  similar  problems. 

This  study  revealed  that  both  ABAQUS  and  ES-DYNA3D  had  difficulties  in  modeling 
the  large  deformations  involved  in  the  gouging  process.  This  was  attributed  to  an 
inherent  limitation  of  the  Fagrangian  mesh  used  by  these  software  packages.  CTH  was 
detennined  to  be  the  best  tool  available  for  the  numerical  investigation  of  hypervelocity 
gouging. 

The  gouging  model  developed  by  Faird  was  further  refined  by  Szmerekovsky 
[14]  to  relate  numerical  results  directly  back  to  experimental  test  conditions  at  the 
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HHSTT.  A  dimensional  analysis  and  similitude  study  was  conducted  to  facilitate  the 
creation  of  CTH  models  based  on  actual  test  sled  configurations.  In  addition  to 
establishing  which  physical  dimensions  and  characteristic  parameters  influence  the 
initiation  of  gouging,  Szmerekovsky  also  considered  factors  for  the  mitigation  of 
gouging,  specifically  the  use  of  coatings.  His  results  reflect  the  successful  use  of  coatings 
at  the  HHSTT  where  an  epoxy  coating  is  applied  to  the  rail  in  the  region  where  the  test 
sled  is  expected  to  reach  its  highest  velocities  as  a  gouging  mitigation  measure. 

Problem  Statement  and  Objectives 

At  the  HHSTT,  a  significant  portion  of  the  rail  is  not  coated  and  observations  of 
rocket  sled  tests  suggest  that  friction  may  have  a  significant  role  in  the  development  of 
conditions  conducive  to  gouging  formation.  Whether  frictional  effects  can  be  accounted 
for  with  a  CTH  model  is  the  concern.  This  issue  is  tied  to  the  method  used  to  define 
material  interfaces,  and  this  research  effort  explores  the  methods  available  in  CTH. 
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II.  Theory  and  Literature  Review 


Chapter  Overview 

This  chapter  provides  background  theory  concerning  friction  and  hydrocodes, 
with  a  focus  on  CTH  and  its  boundary  layer  algorithm.  Also  presented  in  this  chapter  are 
previous  research  efforts  in  the  fields  of  friction  and  high  energy  impact. 

Friction 

Friction  is  a  force  between  two  surfaces,  parallel  to  the  interface,  and  opposing 
relative  motion.  Four  general  mechanisms  are  presented  by  Bayer  [3]  as  the  basis  for 
friction  between  two  solids:  adhesion,  abrasion,  hysteresis,  and  viscous  losses.  Adhesive 
friction  arises  from  the  force  required  to  shear  bonded  junctions  that  form  between  two 
surfaces.  Similarly,  abrasive  friction  is  associated  with  the  force  required  to  deform  the 
surfaces  or  to  cut  a  groove  or  chip.  Hysteresis  is  a  product  of  repeated-cycle 
deformations  and  takes  into  account  that  real  materials  are  not  perfectly  elastic.  Energy 
is  dissipated  in  each  stress  cycle,  with  the  effect  of  either  asperity  deformation  on  a  micro 
scale  or  deformation  of  the  overall  geometry  of  the  solids.  Friction  due  to  viscous  losses 
occurs  when  a  fluid  or  lubricant  is  placed  between  the  two  surfaces.  The  magnitude  of 
the  friction  force  F  can  be  expressed  as  a  sum  of  these  4  components  and  a  corresponding 
expression  for  the  coefficient  of  friction,  p,  can  be  obtained  by  dividing  through  by  the 
nonnal  force,  N: 

^  -  ^ad  +  Fab  +  ^hys+  ^vis  n  It 
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(2.2) 


_  F 

fi  ~  ^  -  ^ad  +  ^ab  +  lihys+  ^vis 

where  subscripts  ad,  ab,  hys,  and  vis  are  the  components  of  friction  relating  to  adhesion, 
abrasion,  hysteresis,  and  viscous  losses  respectively. 

Each  component  of  friction  contributes  in  varying  quantities  that  are  dependent  on 
the  physical  scenario,  but  they  result  in  the  same  overall  effect.  This  effect  is  the 
dissipation  of  energy,  a  most  significant  aspect  of  friction.  The  energy  dissipating  effect 
of  friction  occurs  through  the  generation  of  heat  and  material  deformation.  Bayer  asserts 
that  the  dominant  process  is  the  generation  of  heat,  which  accounts  for  over  90%  of  the 
energy  dissipated.  While  the  4  components  of  friction  are  conceptually  unique,  the 
process  of  identifying  the  specific  contribution  of  each  component  in  a  physical  scenario 
is  often  difficult.  In  a  hypervelocity  gouging  scenario,  it  would  be  reasonable  to  expect 
that  only  the  mechanisms  of  adhesion,  abrasion,  and  viscosity  contribute  to  the  friction 
force. 

Experimental  work  was  conducted  by  Bowden  and  Freitag  [6]  in  1958  to 
characterize  friction  between  metal  surfaces  at  velocities  up  to  800  m/s.  In  these 
experiments,  coefficients  of  friction  were  established  for  steel  on  copper,  aluminum, 
duralumin,  bismuth,  antimony,  and  molybdenum.  The  experiment  made  use  of  a  steel 
sphere  suspended  in  a  magnetic  field  and  accelerated  by  a  rotating  magnetic  field.  The 
freely  spinning  sphere  was  then  decelerated  through  contact  with  three  symmetrically 
arranged  flat  plates  of  the  dissimilar  material  of  interest.  The  friction  force  was 
calculated  by  measuring  the  deceleration  of  the  sphere.  The  friction  coefficient  was  then 
calculated  by  dividing  the  calculated  friction  force  by  the  known  normal  force  applied  by 
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the  three  flat  plates.  While  actual  friction  coefficients  were  highly  dependent  on  the 
materials  involved,  the  results  showed  a  clear  trend  of  decreasing  coefficient  with 
increasing  velocity  in  all  cases.  The  data  fitted  curves  also  suggested  some  minimum 
coefficient  of  friction  greater  than  zero  was  reached  asymptotically  that  varied  from  as 
high  as  .20  for  steel  on  aluminum  to  as  low  as  .05  for  steel  on  bismuth. 

Hydrocodes 

The  primary  tool  used  in  this  research  is  CTH,  a  hydrocode  produced  by  Sandia 
National  Labs.  A  hydrocode  is  a  software  package  developed  to  numerically  simulate 
highly  dynamic  events  involving  shocks.  Anderson  [1]  presents  an  infonnative  overview 
of  this  class  of  computer  programs  which  are  capable  of  simulating  the  propagation  of 
shockwaves  and  computing  properties  such  as  velocities,  stresses,  strains,  and  energy  as 
functions  of  space  and  time.  To  accomplish  this,  hydrocodes  solve  a  discretized  set  of 
the  conservation  equations  of  mass,  momentum,  and  energy,  an  equation  of  state,  a 
material  constitutive  model,  and  a  material  failure  model.  Two  fundamental  approaches 
to  the  solution  of  the  conservation  equations  exist:  the  material  based  Lagrangian 
approach  and  the  spatially  based  Eulerian  approach.  This  historically  resulted  in  two 
distinct  subclasses  of  hydrocodes  whereas  modern  hydrocodes  such  as  CTH  take  a  hybrid 
approach  to  realize  the  benefits  of  both  the  Lagrangian  and  the  Eulerian  methods. 

Lagrangian  hydrocodes  employ  a  computational  grid  that  is  fixed  to  the  material 
with  initial  connections  between  the  grid  points  fonning  cells  of  constant  mass.  The 
computational  grid  need  only  be  generated  where  material  actually  exists  in  the  defined 
problem  space.  The  grid  points  are  constrained  to  move  with  the  local  material  velocity 
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and  values  such  as  velocity,  pressure,  density,  and  temperature  are  calculated  for  each 
grid  point  at  each  time  step  throughout  the  solution  process.  The  immediate  advantage  of 
this  approach  is  the  ability  to  clearly  identify  material  boundaries  and  interfaces.  This 
facilitates  the  process  of  defining  boundary  conditions  at  free  surfaces  as  well  as  contact 
surfaces.  The  frictional  effects  at  the  contact  surfaces  can  also  be  readily  captured. 
Significant  disadvantages  of  a  material  fixed  grid  reveal  themselves  in  large  deformation 
problems  as  time  progresses.  Issues  of  grid  distortion  and  cell  compression  can  produce 
computing  difficulties  as  well  as  significant  errors.  The  extreme  case  of  grid  distortion 
where  negative  area/volume  cells  result  will  halt  the  solution  process.  In  a  hypervelocity 
impact  scenario  this  could  manifest  very  early  (first  few  microseconds)  in  the  solution 
process  due  to  large  defonnations  at  the  contact  surface.  Extreme  cell  compression 
forces  the  hydrocode  to  reduce  the  time  step  resulting  in  an  overall  computational  time 
that  may  not  be  acceptable  or  feasible.  To  proceed  with  the  solution  would  require  a  new 
computational  grid  to  be  overlaid  with  properties  calculated  from  the  distorted  grid, 
keeping  in  mind  the  need  to  conserve  mass,  momentum,  and  energy.  Material  mixing  in 
the  impact  region  would  also  not  be  captured  by  a  Lagrangian  hydrocode  since  mixed 
material  cells  are  not  possible. 

Eulerian  hydrocodes  avoid  the  problems  of  grid  distortion  and  cell  compression 
through  the  use  of  a  spatially  fixed  grid  which  results  in  fixed  volume  cells.  However, 
the  computational  grid  must  be  defined  for  the  entire  problem  space,  including  originally 
empty  regions  to  which  material  may  propagate  as  the  solution  progresses.  This 
approach  computes  the  flow  of  mass,  momentum,  and  energy  across  the  cell  boundaries. 
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Since  mass  flow  through  the  cells  is  permitted,  this  approach  accommodates  problems 
involving  material  diffusion  and  mixing.  The  mixing  of  material  at  contact  surfaces  can 
be  accurately  modeled.  However,  this  ability  to  handle  mixed  cells  can  result  in 
ambiguous  material  boundaries  and  interfaces.  No  distinct  material  boundaries  exist 
within  a  given  cell,  complicating  the  task  of  applying  boundary  conditions.  To  reduce  the 
impact  of  this,  smaller  grid  point  spacing  can  be  applied  at  the  expense  of  increased 
computation  time.  With  no  distinct  material  interfaces,  the  study  of  friction  on  contact 
surfaces  becomes  more  arduous  and  less  accurate. 

CTH 

The  hydrocode  CTH  was  developed  at  Sandia  National  Laboratories  to  handle  3- 
D,  multi-material,  large  deformation,  strong  shock  problems  [11].  CTH  employs  an 
Eulerian  mesh  where  all  quantities  are  cell  centered  except  for  the  face  centered 
velocities.  The  two-step  Eulerian  solution  scheme  used  by  CTH  involves  a  Lagrangian 
step  followed  by  a  remap  step  and  is  second  order  accurate. 

In  the  Lagrangian  step,  the  mesh  is  allowed  to  distort  with  the  material.  During 
this  step,  conservation  of  volume,  mass,  momentum,  and  energy  are  conserved.  Since  the 
mesh  is  fixed  to  the  material  in  this  step,  mass  is  conserved  automatically.  The  governing 
equations  are  replaced  with  explicit  finite  volume  approximations.  However,  solution  of 
these  equations  requires  a  timestep  small  enough  to  limit  the  volume  change  in  individual 
cells  and  to  limit  the  distance  a  wave  can  travel  to  less  than  one  cell  width. 

The  second  step  remaps  the  distorted  cells  back  to  the  original  Eulerian  mesh. 

The  volume  flux  between  cells  is  calculated  first  based  on  the  cell  face  velocities.  A 
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high-resolution  interface  tracker  algorithm  is  then  implemented  to  ensure  the  right 
materials  are  moved  with  the  volume.  Each  material’s  corresponding  mass,  momentum, 
and  energy  are  then  moved. 

CTH  relies  heavily  on  graphics  to  present  results  to  the  user  but  data  from  tracer 
points  can  also  be  obtained.  These  tracer  points,  which  can  be  Eulerian  or  Lagrangian, 
need  to  be  defined  before  the  problem  is  submitted  to  CTH  for  solution.  Some 
knowledge  of  the  expected  regions  of  interest/activity  is  necessary  for  the  strategic 
placement  of  these  points.  Data  for  the  rest  of  the  defined  mesh  can  be  obtained  through 
the  use  of  contour  and  vector  plots.  The  data  range  for  such  plots  is  defined  by  the  user, 
so  some  knowledge  of  the  expected  results  along  with  trial  and  error  is  necessary  to 
generate  meaningful  results  with  the  graphical  post  processing  capability. 

Boundary  Layer  Algorithm 

CTH  was  originally  developed  to  simulate  problems  involving  strong  shock 
waves  where  pressures  are  very  large  but  shear  stresses  are  relatively  small.  Under  this 
premise,  two  simplified  options  were  implemented  for  the  treatment  of  effects  at  material 
interfaces.  The  default  option  is  to  assume  materials  are  joined  at  their  contact  surface  by 
requiring  the  velocity  be  constant  across  the  interface.  The  second  option  is  the  slide  line 
approach.  Before  further  discussion  regarding  this  approach,  it  should  be  noted  that  the 
state  of  stress  at  a  point  can  be  represented  by  the  Cauchy  stress  tensor  which  can  be 
expressed  in  terms  of  a  spherical  component  and  a  deviatoric  component  as  follows  [12]: 

(2.3) 


°"ll  °"l2  °"l3  °"m  0  0  °"ll  °"l2  °"l 

&12  ®22  ®22  =  0  <Jm  0  +  CT12  O' 22  CT  2 

CT  i  CT  o  i  CT 11  0  0  CT  CT 1 1  CT  ^  ->  CT  ■> 
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where  crm  is  the  mean  stress  or  pressure  and  is  define  as: 


cr  = 


(Tn  "I"  Cr22  ^33 


(2.4) 


Continuing  on,  the  alternate  slide  line  approach  designates  the  contact  surface  as  a 
frictionless  sliding  interface  by  setting  the  material  strengths,  and  hence  the  deviatoric 
stress,  of  mixed  material  cells  to  zero.  Only  pressure  is  supported  by  these  mixed  cells. 
Even  with  a  model  that  initially  aligns  the  material  surfaces  with  the  grid  lines,  as  the 
materials  make  contact  and  deform,  the  material  interfaces  will  shift  into  the  cell  interiors 
and  result  in  mixed  material  cells  as  shown  in  Figure  2.1.  The  slide  line  approach  can 
have  the  effect  of  reducing  the  strength  of  a  projectile  surface  since  that  surface  is  treated 
as  a  fluid.  This  fluid  layer  can  be  swept  away  as  the  projectile  penetrates  a  target 
resulting  in  a  distorted  projectile  shape.  This  was  an  acceptable  result  for  the  impact 
studies  of  interest  when  CTH  was  originally  developed,  since  erosion  of  the  projectile 
surface  and  deformation  of  the  projectile  in  general  as  it  penetrated  the  target  was 
experimentally  observed.  However,  interest  in  lower  velocity  penetration  and  perforation 
problems,  below  1.5  km/s,  have  become  a  subject  of  interest  and  in  these  scenarios, 
significant  deformation  of  the  projectile  is  not  physically  observed.  CTH  simulations  of 
these  low  velocity  penetration  events  using  the  slide  line  approach  can  produce  erroneous 
results  by  showing  projectile  defonnation  which  may  not  reflect  physical  phenomena. 

In  an  attempt  to  more  accurately  simulate  low  speed  penetration  and  perforation 
problems,  a  boundary  layer  algorithm  was  developed  to  improve  the  way  CTH  handles 
sliding  interfaces  [13].  Since  this  algorithm  was  developed  for  these  particular  problems, 
its  implementation  in  CTH  is  limited  in  2D  to  axisymmetric  problems.  The  basic 
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approach  of  the  boundary  layer  algorithm  is  to  move  the  sliding  interface  into  the  target 
material.  This  eliminates  the  purely  numerical  effect  of  erosion  of  the  penetrator  surface. 

The  first  phase  of  the  boundary  layer  algorithm  identifies  cells  in  two  distinct 
layers,  a  hard  boundary  layer  and  a  soft  boundary  layer.  Within  these  boundary  layers 
lies  the  interface  layer  which  contains  the  actual  contact  surface.  Within  the  soft 
boundary  layer  lies  the  slip  layer  which  is  the  artificial  sliding  surface. 

The  interface  layer  is  defined  first.  CTH  calculates  the  material  volume  fractions, 
(f)h  and  (j)5 ,  corresponding  to  the  hard  and  soft  materials  respectively,  for  each  mesh  cell. 

To  locate  the  material  interface,  the  components  of  the  vector  gradients  V<f>h  and  Vtf)s  are 
then  calculated  using  central  difference  formulas: 


2Ax 


»  (VAy)iJ  = 


(0h)i+l,j 

2Ay 


(2.5) 


wju  = 


(<f>s  )MJ  /V7^  ,  -(^),  I., 


2Ax 


.  wju  = 


2A y 


(2.6) 


in  which  Ax  and  Ay  are  the  cell  widths  and  the  subscripts  i,j  indicate  the  cell  coordinates. 
Figure  2.1  illustrates  a  typical  mesh  with  cell  coordinates  for  a  representative  mesh  region 
indicated.  Cells  in  which  both  the  |(V^/i);.  y.|  >  0.1  and  |(V$s)iy.|  >  0.1  define  the  interface 


layer.  This  normally  results  in  a  layer  that  is  two  to  three  cells  thick  depending  where  the 
physical  material  interface  is  in  relation  to  the  cell  edges.  Figures  2.2  and  2.3  illustrate 
the  cases  where  the  physical  material  interface  does  or  does  not  line  up  with  cell  edges. 
The  latter  case  establishes  mixed  material  cells  from  the  start.  An  effective  unit  normal 
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vector  ny,  which  is  perpendicular  to  the  material  interface  and  points  into  the  hard 
material,  is  also  calculated  for  each  cell  in  the  interface  layer  as  shown  in  Figures  2.2  and 
2.3: 


<y*k)u 

(v«U,y -WXj 


(2.7) 


x ,  1 


Figure  2.1  Typical  Mesh  Showing  Sample  Cell  Coordinates 
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Figure  2.2  Interface  Layer,  Material  Interface  Aligned  With  Cell  Edge 
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Figure  2.3  Interface  Layer,  Material  Interface  Inside  Cells 


The  two  boundary  layers  are  defined  next.  A  cell  is  in  one  of  the  boundary  layers 
if  its  cell  center  is  within  some  user  defined  distance,  ww,  of  an  interface  layer  cell.  If  the 
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cell  also  has  a  soft  material  volume  fraction,  <f>s ,  of  at  least  .99,  than  it  is  in  the  soft 

boundary  layer.  Otherwise,  it  is  in  the  hard  boundary  layer.  This  effectively  designates 
almost  all  mixed  material  cells  as  hard  boundary  layer  cells. 

Finally,  the  slip  layer  is  defined  as  those  cells  in  the  soft  boundary  layer  that  are 
also  within  some  user  defined  distance,  wsi,  of  an  interface  layer  cell.  The  slip  layer 
models  a  frictionless  sliding  interface  by  setting  the  deviatoric  stress  to  zero. 

Shown  in  figures  2.4  and  2.5  are  the  boundary  layer  configurations  for  the  cases 
where  the  physical  material  interface  does  or  does  not  align  with  the  cell  edges.  In  the 
latter  case,  the  hard  boundary  layer  actually  contains  the  mixed  material  cells,  effectively 
increasing  the  size  of  the  penetrator.  These  mixed  cells  also  have  their  yield  stress  set  to 
the  hard  material  yield  stress  to  prevent  nonphysical  erosion  of  the  hard  material  as 
happens  with  the  slide  line  method. 


Figure  2.4  Boundary  Layers,  Material  Interface  Aligned  With  Cell  Edge 
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The  second  phase  of  the  Boundary  Layer  Algorithm  calculates  friction  forces  for 
cells  that  exist  in  either  Boundary  Layer  but  not  the  Interface  Layer,  which,  in  general, 
contain  mixed  material  cells.  Friction  forces  are  not  applied  to  Interface  Layer  cells  due 
to  limitations  in  CTH’s  ability  to  handle  mixed  material  cells.  Friction  is  not  applied  as  a 
direct  force  which  would  contribute  to  the  deviatoric  stress  but  rather  as  a  body  force  that 
contributes  to  the  momentum  balance.  This  friction  body  force  is  computed  as  follows. 

For  each  cell  in  one  of  the  boundary  layers,  an  average  velocity  of  the  hard 
material,  Vh,  and  an  average  velocity  of  the  soft  material,  Vs,  is  calculated.  Contributing 
values  are  typically  taken  from  cells  within  an  approximate  5  cell  radius,  disregarding 
slip  layer  cells  which  may  not  accurately  reflect  the  bulk  motion  of  the  soft  material.  The 
relative  velocity  vector  of  hard  and  soft  material  for  each  cell  is  then  evaluated  from: 

V  =  Vh-Vs  (2.8) 
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Similarly  the  average  unit  normal  vector,  n  ,  and  the  average  Cauchy  stress  tensor,  T  ,  for 

each  cell  are  calculated  from  the  same  contributing  cells.  Since  V  may  not  necessarily 

be  tangential  to  the  contact  surface,  a  unit  vector,  u,  which  is  tangential  to  the  contact 

surface  and  points  in  the  direction  of  the  friction  force  acting  on  the  soft  material,  as 

shown  in  Figure  2.6,  is  calculated: 

(nxV)xn 
(nxV)xn 


Figure  2.6  Diagram  of  Unit  Vector  u 


The  shear  traction  vector  t/is  then  calculated  from 

tf=  +/min{0  ,  n -(Tn)  }u  (2.10) 

where/is  the  user  defined  coefficient  of  friction  and  the  scalar  quantity  n  ■  (Tn)  is  the 
component  of  stress  in  the  direction  of  n  ,  also  referred  to  as  the  normal  traction,  and  can 
also  be  expressed  in  matrix  notation  as: 

CTn  cr12  cr13  lx 

n-(Tn)  =  %  l2  /”]•  an  a12  a23  /”  (2.11) 

<TU  <J23  (733  l3 

For  the  shear  traction  to  be  nonzero,  the  normal  traction  must  be  compressive  (negative 
value).  The  sign  is  negative  if  the  cell  is  in  the  soft  layer  to  produce  a  positive  shear 
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traction  acting  on  the  soft  material  in  the  direction  of  u.  The  sign  is  positive  if  the  cell  is 
in  the  hard  layer  to  produce  a  negative  shear  traction  acting  on  the  hard  material 
opposition  the  direction  of  u.  The  shear  traction  is  then  converted  to  a  body  force 
density: 

b/=  t//  wbi  (2.12) 

which  appears  in  the  momentum  balance  equation: 

pa  =  V-T  +  b/  (2.13) 

in  which  p  is  the  density,  a  is  the  acceleration,  and  T  is  the  Cauchy  stress  tensor. 

Due  to  issues  of  mesh  coarseness  and  uncertainty  in  the  contact  surface  shape 
which  can  affect  the  number  of  cells  included  in  the  interface  layer,  the  amount  of  friction 
calculated  at  any  given  will  fluctuate.  Over  time  however,  this  fluctuation  is  averaged 
out  and  the  net  friction  force  applied  to  the  solution  approaches  the  correct  amount. 

To  apply  the  boundary  layer  algorithm,  the  user  must  designate  a  hard  material 
and  a  soft  material  for  each  pairing  of  materials  that  may  come  into  contact  and  require  a 
boundary  layer  interface  to  be  defined.  The  width  of  the  boundary  layers  can  be  defined 
by  the  user  but  is  not  necessary.  Also,  to  explicitly  add  friction  to  a  scenario,  a 
coefficient  of  friction  must  be  defined  by  the  user. 

An  option  also  is  available  to  artificially  increase  the  yield  stress  of  the  hard 
material  to  prevent  non  physical  plastic  defonnation  of  a  projectile.  This  can  occur  if  the 
applied  shear  stress  is  close  enough  to  the  material  yield  stress  for  numerical  noise  to  be  a 
concern. 


19 


The  boundary  layer  algorithm  also  restricts  material  in  the  slip  layer  from  flowing 
counter  to  the  direction  of  the  projectile  which  the  user  can  override.  The  velocity 
component  of  material  in  the  Slip  Layer  parallel  to  the  projectile  velocity  vector  is  not 
pennitted  to  be  less  than  the  projectile  velocity.  If  counterflow  is  permitted  in  the  slip 
layer,  the  algorithm  can  fail  to  compute  that  the  projectile  ever  comes  to  rest  since  target 
material  would  continually  be  ejected  from  the  slip  layer. 

To  validate  the  improvements  made  to  CTH  with  the  implementation  of  the 
boundary  layer  algorithm.  Silling  [13]  developed  several  test  cases  to  compare  CTH 
solutions  to  available  experimental  data.  One  case  of  interest  simulates  the  deep 
penetration  of  an  aluminum  target  by  a  hemispherical-nosed  maraging  steel  rod.  This 
case  was  experimentally  investigated  by  Forrestal,  Brar,  and  Luk  [7].  The  projectile  was  a 
T200  maraging  steel  rod  of  length  7 1 . 12mm  and  diameter  7. 1 1mm.  The  target  was  a 
6061-T651  aluminum  cylinder  of  length  198mm  and  diameter  152mm.  Material 
properties  are  presented  in  Table  2.1,  where  Y  is  the  tensile  yield  stress,  E  is  the  modulus 
of  elasticity,  and  v  is  Poisson’s  ratio.  The  impact  velocity  was  measured  at  959m/s.  For 
the  CTH  simulations,  Silling  employed  a  2D  axisymmetric  model  that  contained  the 
entire  target  and  projectile.  The  projectile  was  modeled  as  an  elastic  perfectly  plastic 
material  but  the  target  was  modeled  with  the  Johnson-Cook  viscoplastic  model.  Both 
materials  were  modeled  using  a  Mie-Gruneisen  equation  of  state. 
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Table  2. 1  Material  Properties,  Forrestal,  Brar,  and  Luk  Experiment 


T200  Maraging  Steel 

6061-T651  Aluminum 

Y 

1400  MPa 

240  MPa 

E 

200  Gpa 

68  GPa 

V 

.31 

.33 

Experimental  results  indicated  a  penetration  depth  of  109mm.  Inspection  of  the 
target  material  also  revealed  a  thin  layer  along  the  inside  of  the  resulting  cavity  that 
suggested  the  material  had  melted.  The  formation  of  the  melted  layer  was  assumed  to  be 
a  frictional  effect  resulting  from  the  contact  between  the  target  and  the  nose  of  the 
projectile.  Silling  obtained  solutions  using  mesh  cell  sizes  of  .08cm  and  .04cm  and  a 
friction  coefficient  of  .06.  The  .08cm  mesh  under  predicted  the  depth  at  102mm  while 
the  .04cm  mesh  over  predicted  the  depth  at  1 15mm. 
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ITT.  Methodology 


Chapter  Overview 

Presented  in  this  chapter  are  the  methods  used  to  investigate  the  three  methods 
employed  by  CTH  to  model  material  interfaces.  The  three  methods  will  be  referred  to  as 
1)  No  Slide  Line,  2)  Slide  Line,  and  3)  Boundary  Layer.  The  first  method  is  the  default 
and  assumes  materials  are  joined  at  their  contact  surfaces.  The  Slide  Line  method 
establishes  a  frictionless  sliding  interface  at  the  contact  surfaces.  The  Boundary  Layer 
method  implements  the  boundary  layer  algorithm  where  the  sliding  interface  is 
established  in  the  target  material.  These  three  methods  are  compared  in  a  nonnal  impact 
scenario  and  a  tangential  sliding  scenario. 

Normal  Impact,  Part  1 

The  deep  penetration  model  developed  by  Silling  [13]  to  validate  the  boundary 
layer  algorithm  is  used  for  the  initial  phase  of  this  research  and  is  illustrated  in  Figure  3.1. 
This  axisymmetric  model  was  defined  in  2D  cylindrical  coordinates  with  a  symmetry 
boundary  condition  at  x  =  0  and  sound  speed  absorbing  boundary  conditions  at  the  three 
remaining  mesh  boundaries.  The  sound  speed  absorbing  boundary  condition  simulates  a 
semi  infinite  boundary  to  prevent  shock  reflections  back  into  the  mesh  [4],  This  model 
was  used  to  investigate  four  cases: 

1 .  The  model  as  provided,  which  employs  the  Boundary  Layer  method,  to  verify 
reproducibility  of  the  reported  results 

2.  The  model  employing  the  Slide  Line  method 

3.  The  model  employing  the  No  Slide  Line  method 
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4.  The  Boundary  Layer  model  with  friction  coefficient  (p)  changed  from  0.02  to  0.0 


Figure  3.1  Graphic  of  Impact  Scenario 


These  impact  cases  are  studied  to  gain  insight  on  how  CTH  handles  material 
interfaces.  The  results  are  used  to  quantify  the  improvements  that  the  boundary  layer 
algorithm  provides.  The  fourth  case  is  expected  to  isolate  the  effect  of  shifting  the  sliding 
interface  away  from  the  actual  material  interface  while  still  assuming  frictionless  sliding. 
These  solutions  were  resolved  out  to  300  microseconds  to  ensure  maximum  penetration 
depths  were  achieved  for  comparison. 
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Lagrangian  tracers  are  used  to  capture  the  conditions  which  develop  at  specific 
points  along  the  contact  surface  between  the  projectile  and  target.  The  use  of  Lagrangian 
tracers  is  necessary  since  the  contact  surface  moves  with  the  projectile. 

The  original  .04cm  mesh  was  also  further  refined  to  a  mesh  cell  size  of  ,01cm. 
This  was  accomplished  due  to  numerical  artifacts  observed  at  the  original  .04cm  cell  size. 
Solutions  were  obtained  again  using  all  three  methods  for  defining  material  interfaces. 
The  fourth  case  with  coefficient  of  friction  set  to  0.0  was  not  solved  with  the  refined 
mesh.  The  two  Boundary  Layer  cases  were  sufficiently  similar  with  the  ,04cm  mesh  that 
further  investigation  was  not  warranted.  Also,  since  the  boundary  layer  algorithm  does 
not  support  data  parallel  computing,  solutions  were  only  resolved  out  to  50  microseconds 
since  the  finer  meshes  were  significantly  more  computationally  intensive  and  time 
consuming  running  on  a  single  workstation.  These  solutions  were  also  generating  output 
data  files  approaching  the  2GB  size  limitation  imposed  by  CTH. 

Normal  Impact,  Part  2 

The  deep  penetration  model  was  then  modified  to  incorporate  actual  materials 
used  at  the  HHSTT.  The  penetrator  material  was  changed  to  a  strain  rate  independent 
Steinberg-Guinan-Lund  viscoplastic  model  of  Vascomax  300.  The  target  material  was 
changed  to  a  Johnson-Cook  viscoplastic  model  of  Iron,  the  most  similar  material  model 
to  1080  steel  available.  Tabular  equation  of  state  data  was  also  available  for  both 
materials.  The  material  properties  are  presented  in  Table  3.1,  where  Tmeit  is  the  melting 
temperature.  Solutions  for  all  three  material  interface  methods  were  obtained  with  a 
mesh  size  of  .01cm.  Since  the  penetration  depths  encountered  with  the  change  in 
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materials  was  significantly  less,  the  overall  mesh  dimensions  were  reduced  to  allow 
solution  times  out  to  100  microseconds. 


Table  3.1  Material  Properties,  HHSTT 


Vascomax  300 

1080  Steel 

Iron 

Y 

2000  MPa 

700  MPa 

572  MPa 

E 

200  GPa 

200  GPa 

172  GPa 

V 

.31 

.25 

.28 

Tmelt 

1650  K 

1835  K 

1811  K 

Pure  Sliding 

The  final  phase  of  this  study  attempts  to  isolate  the  frictional  effects  of  a  high 
energy  impact  by  simulating  a  pure  sliding  scenario.  The  normal  impact  model  used 
previously  is  modified  by  removing  the  material  directly  in  the  path  of  the  projectile. 
This  models  a  cylinder  sliding  inside  a  cylinder  and  is  illustrated  in  Figure  3.2.  The  only 
contact  surface  in  this  scenario  is  parallel  to  the  direction  of  motion  however  there  is 
continuous  contact  along  the  length  of  the  projectile.  Shear  forces  on  that  surface  should 
be  purely  a  result  of  sliding  friction.  Solutions  for  a  projectile  velocity  of  lkm/s  were 
resolved  out  to  50  microseconds. 
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Figure  3.2  Graphic  of  Sliding  Scenario 
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IV.  Analysis  and  Results 


Normal  Impact  Results,  Part  1 

The  initial  solutions  using  the  normal  impact  model,  with  initial  velocity  of 
959m/s  and  mesh  cell  size  of  .01cm,  support  the  assertion  that  the  Boundary  Layer 
method  produces  more  realistic  solutions  in  the  simulation  of  deep  penetration.  The 
maximum  penetration  depths  for  the  four  cases  studied  are  presented  in  Table  4.1,  along 
with  the  experimental  result  reported  by  Forrestal,  et  al  [7]  and  the  CTH  result  reported 
by  Silling  [13].  It  should  be  noted  that  the  maximum  penetration  depths  are  reached  prior 
to  300  microseconds  and  the  final  depth  computed  is  marginally  less.  The  discrepancy  in 
the  maximum  depth  and  the  final  position  could  be  a  result  of  elastic  recovery  in  the 
materials  or  a  numerical  artifact.  The  exact  penetration  depth  reported  by  Silling  could 
not  be  repeated  but  this  is  most  likely  a  result  of  differences  in  the  version  of  CTH  that 
was  used.  The  current  version  of  CTH  is  6.01  while  the  boundary  layer  algorithm  was 
implemented  in  version  1.027  although  there  is  no  indication  of  changes  made  to  the 
boundary  layer  algorithm  since  its  implementation.  The  significance  of  these  results  is 
that  they  support  the  assertion  that  the  Boundary  Layer  method  offers  much  improvement 
over  the  other  two. 


Table  4. 1  Penetration  Depths  (.04cm  ’ 

Mesh) 

Experiment 

Silling 

Boundary 

Layer 

p  =  .06 

Slide 

Line 

No  Slide 
Line 

Boundary 

Layer 

p  =  0. 

Penetration  Depth  (mm) 

109 

115 

99.1 

78.2 

80.5 

103.6 

%  Difference  from 

Experiment 

+6% 

-9% 

-28% 

-26% 

-5% 
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Of  interest  also  are  the  numerous  “thermodynamic  warnings”  that  are  generated 
with  the  Slide  Line  method.  These  warnings  are  generated  when  the  equation  of  state 
model  calculates  an  unrealistic  thermodynamic  state  for  any  cell,  which  usually  manifests 
as  a  negative  temperature  value.  The  user  sets  an  arbitrary  limit  for  how  many  of  these 
warnings  can  occur  before  the  CTH  calculation  aborts.  While  there  is  no  limit  as  to  how 
high  this  limit  can  be,  the  accuracy  of  the  solutions  may  be  questionable  if  a  significant 
number  of  these  thermodynamic  warnings  are  generated. 

Pressure  contour  plots  for  these  4  cases  are  presented  in  Figures  4.1  -  4.4  with 
black  dots  indicating  locations  of  material  tracers.  These  tracers  provide  a  visual 
indication  of  local  material  movement.  At  50  microseconds,  the  impact  events  have 
developed  enough  to  distinguish  unique  characteristics  of  each  approach.  Pressure  bands 
of  alternating  intensity  are  visible  in  the  projectile  while  a  relatively  smooth  decreasing 
pressure  gradient  exists  in  the  target.  This  is  suspected  to  be  a  result  of  the  dimensions  of 
the  model.  The  relatively  small  projectile  dimensions  allow  pressure  waves  to  reflect  at 
the  free  surfaces  and  interact.  The  target  however  is  large  enough  that  pressure  waves  are 
not  observed  to  reflect  off  the  side  and  bottom  free  surfaces.  The  band  of  negative 
pressure  near  the  target’s  top  free  surface  does  indicate  some  boundary  effect  but  was 
limited  to  the  immediate  vicinity  of  the  surface. 
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Figure  4.1  Pressure  Contour  at  50  jjs,  No  Slide  Line 
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Figure  4.2  Pressure  Contour  at  50  gs,  Slide  Line 
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Figure  4.3  Pressure  Contour  at  50  gs,  Boundary  Layer,  g=0.06 
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Figure  4.4  Pressure  Contour  at  50  gs,  for  Boundary  Layer,  g=0.0 
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Additionally,  these  figures  clearly  show  the  extent  to  which  preservation  of  the 
projectile  shape  is  maintained  through  use  of  the  Boundary  Layer  method.  However,  one 
disconcerting  feature  which  is  most  prevalent  with  the  Boundary  Layer  method  is  the 
existence  of  horizontal  material  bands  between  the  projectile  and  the  target.  This  can  be 
attributed  to  the  coarseness  of  the  mesh  which  makes  it  more  difficult  for  CTH  to  track 
material  surfaces.  Small  distortions  of  the  projectile  surface  can  affect  which  cells  are 
identified  as  boundary  layer  cells  as  can  be  seen  in  Figure  4.5  with  the  boundary  layer 
type  values  defined  in  Table  4.2.  Since  the  target  material  in  mixed  cells  is  attached  to 
the  projectile,  additional  target  material  is  pulled  along  with  the  projectile  whenever  the 
hard  boundary  layer  shifts  towards  the  target  material.  The  effect  of  the  frozen  slip  layer 
compounds  this  issue  since  that  material  is  forced  to  move  down  into  the  target  with  the 
projectile.  This  effect  is  expected  to  be  reduced  as  mesh  size  is  reduced  but  not  entirely 
eliminated.  (It  should  be  noted  that  the  program  used  to  generate  these  plots,  namely 
CTHPLT,  applies  a  smoothing  technique  that  can  produce  misleading  results.  Since  these 
are  contour  plots,  the  program  does  not  allow  discreet  jumps  between  the  data  values;  all 
plotted  data  must  transition  smoothly.  In  Figure  4.5,  no  soft  boundary  layer,  shown  in 
red,  actually  exists  in  the  hard  projectile  material;  this  is  purely  a  result  of  smoothing.) 
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Figure  4.5  Boundary  Layer  Type  at  50  ps,  p=0.06 


AMSEOS 


Table  4.2  Boundary  Layer  Types 


0 

Cell  is  not  in  a  boundary  layer 

1 

Cell  is  in  a  soft  boundary  layer  but  not  a  slip  layer 

2 

Cell  is  in  a  hard  boundary  layer 

3 

Cell  is  in  a  non-frozen  slip  layer 

4 

Cell  is  in  a  frozen  slip  layer 

The  temperature  profiles  are  also  useful  in  discussing  the  differences  in  these 
solutions  and  are  shown  in  Figures  4.6  -  4.9.  The  No  Slide  Line  and  the  Slide  Line  cases 
show  similar  elevated  temperature  profiles  while  the  Boundary  Layer  cases  show 
relatively  little  change  in  temperature.  An  increase  in  temperature  can  be  interpreted  as 
localized  plastic  deformation,  since  these  solutions  do  not  account  for  heat  transfer. 
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Considering  the  significant  deformation  resulting  from  the  impact  event,  elevated 
temperatures  near  the  impact  surface  would  be  expected. 

It  should  be  noted  at  this  point  that  these  observations  and  assertions  are  based  on 
a  solution  obtained  with  a  relatively  coarse  mesh.  A  finer  mesh  is  necessary  before  any 
realistic  comparisons  between  the  material  interface  models. 
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Figure  4.6  Temperature  at  50  ps,  No  Slide  Line 
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Figure  4.7  Temperature  at  50  ps,  Slide  Line 
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Figure  4.8  Temperature  at  50  ps,  Boundary  Layer,  p=.06 
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Figure  4.9  Temperature  at  50  ps,  Boundary  Layer,  g=0.0 


The  second  set  of  solutions  obtained  using  a  mesh  cell  size  of  .01cm  appears  to 
negate  the  assertions  that  were  previously  made  with  a  mesh  cell  size  of  .04cm.  Figures 
4.10  -  4.12  again  show  the  pressure  contours  for  the  three  cases  at  50  microseconds. 
Similar  pressure  bands  in  the  projectile  are  present  in  all  three  cases.  Perhaps  the  key 
observation  made  from  these  figures  is  that  the  shape  of  the  penetrator  does  not  distort 
noticeably  in  any  of  the  three  cases.  The  shape  of  the  deformed  target  varies  somewhat 
but  the  differences  seen  along  the  vertical  face  is  most  likely  due  to  mesh  coarseness 
affecting  the  material  interface  tracker  algorithm. 
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Figure  4.10  Pressure  Contour  at  50  gs,  No  Slide  Line 
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Figure  4.12  Pressure  Contour  at  50  ps,  Boundary  Layer,  |j=0.06 


A  projectile  material  tracer  located  in  the  nose  of  the  projectile,  with  initial 
coordinates  x=.205  and  y=.305,  was  selected  to  compare  pressure,  deviatoric  shear  stress, 
and  normal  stress  histories  for  each  method.  The  conditions  at  this  tracer  are 
representative  of  the  overall  projectile  contact  surface  since  its  position  relative  to  the 
projectile  and  the  contact  surface  remains  essentially  constant,  less  than  .002cm 
fluctuation.  These  tracer  history  plots  are  shown  in  Figure  4. 13  -  4.15  and  are  virtually 
indistinguishable.  Target  material  tracers  were  observed  to  displace  beyond  the  right 
edge  of  the  projectile  within  the  first  10  microseconds  due  to  severe  plastic  deformation 
of  the  target  material  along  the  contact  surface.  Since  this  places  them  outside  the  impact 
region,  their  data  history  would  not  reflect  conditions  which  develop  near  the  material 
interface  for  the  entirety  of  the  solution  time. 
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The  penetration  depths  at  50  microseconds  are  compared  in  Table  4.2.  While  the 
depth  varies  by  up  to  .377cm  with  the  .04cm  mesh,  the  max  difference  is  only  .056cm 
with  the  .01cm  mesh. 


Table  4.2  Penetration  Depths  at  50  ps 


Boundary  Layer 

p  =  .06 

Slide  Line 

No  Slide  Line 

.04cm  mesh 

4.225  cm 

3.848  cm 

3.890  cm 

.01cm  mesh 

4.331  cm 

4.367  cm 

4.31 1  cm 

38 


PRESSURE  (I0SD/CM2) 


L  POINT  161 


TIME  (10"6S) 

HAL32:  Steel  ->  Al,  959  m/s.  no  Slide  Line,  dx=0.01 


L  POINT  161 


TIME  (10"6S) 

HAL32:  Steel  ->  Al,  959  m/s.  Slide,  dx=0.01 


a)  No  Slide  Line 


b)  Slide  Line 


L  POINT  161 


TIME  (10"6S) 
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c)  Boundary  Layer 

Figure  4.13  Projectile  Surface  Pressure  History 
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c)  Boundary  Layer 

Figure  4.14  Projectile  Surface  Deviatoric  Shear  Stress  History 
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Figure  4.15  Projectile  Surface  Normal  Stress  History 
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These  results  suggest  that  as  the  size  of  mesh  cells  decrease,  the  differences  in 
solutions  obtained  with  the  three  methods  for  defining  material  interfaces  become 
negligible.  Using  the  Boundary  Layer  method  to  preserve  the  shape  of  the  penetrator 
does  not  appear  to  be  necessary.  At  this  level  of  mesh  refinement,  the  effect  of  numerical 
erosion  of  the  penetrator  surface  is  not  a  significant  factor.  The  only  apparent  advantage 
to  using  the  Boundary  Layer  method  is  the  capability  to  explicitly  add  friction  to  the 
solution  except  that  in  the  current  scenario,  where  our  solution  is  underestimating  the 
maximum  penetration  depth,  adding  friction  results  in  a  solution  even  further  from  the 
experimentally  measured  depth.  This  method  doesn’t  even  appear  to  be  advantageous  for 
the  problem  of  penetration  and  perforation  for  which  it  was  developed.  However,  if  one 
considers  that  the  boundary  layer  algorithm  was  developed  in  1992  and  that  the  deep 
penetration  model  used  to  validate  it  was  only  refined  to  a  .04  cm  mesh,  the  need  for  this 
method  becomes  almost  apparent.  Given  the  level  of  computing  capability  present  in 
1992,  the  computational  requirements  to  obtain  solutions  for  a  .01cm  mesh  versus  a 
.04cm  mesh  were  most  likely  time  and  cost  prohibitive.  Even  with  the  modem  computers 
available  in  the  AFIT/ENY  lab,  a  full  solution  out  to  300  microseconds  was  not  feasible. 
The  benefits  of  the  Boundary  Layer  method  then  would  appear  to  only  apply  for 
relatively  coarse  meshes.  However,  it  would  be  premature  to  discard  this  technique  as  a 
viable  option  based  on  the  limited  results  obtained  to  this  point. 

A  comparison  of  the  temperature  profiles  shown  in  Figures  4.16-4.18  indicates  a 
difference  in  the  solutions  that  may  be  significant.  While  the  total  area  of  increased 
temperature  is  comparable,  the  temperature  magnitudes  are  not.  The  magnitude  of  the 
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temperature  increase  and  associated  plasticity  is  lower  for  the  Boundary  Layer  method  as 
observed  previously  for  the  .04  cm  mesh  solutions.  Also,  since  the  projectile  shows 
virtually  no  deformation,  the  plasticity  is  isolated  to  the  target  material.  Considering  the 
differences  in  material  properties  for  the  steel  projectile  and  the  aluminum  target,  this  is 
not  a  surprising  result. 

These  test  cases  have  served  to  provide  some  insight  concerning  the  effects  of  the 
three  methods  for  defining  material  interfaces  but  a  scenario  with  more  material 
interaction  is  necessary  to  obtain  any  more  useful  observations.  Since  previous  gouging 
simulations  using  the  real  materials  have  shown  that  significant  material  interaction 
occurs  under  the  right  conditions,  the  next  phase  of  this  research  examines  an  impact 
scenario  with  real  materials. 
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Figure  4.16  Temperature  at  50  ps,  No  Slide  Line 
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Before  proceeding  however,  the  issue  alluded  to  previously  of  the  horizontal 
material  bands  between  the  vertical  surfaces  of  the  projectile  and  the  target  was  given 


more  consideration.  At  a  mesh  cell  size  of  .0 1cm,  these  material  bands  are  still  present  in 
the  Boundary  Layer  case  as  can  be  seen  in  the  pressure  and  temperature  contours  above. 
The  assertion  made  previously  was  that  the  fonnation  of  these  bands  is  a  product  of  mesh 
coarseness  and  the  adverse  effect  it  has  on  accurately  defining  boundary  layers.  The 
amount  of  material  banding  was  reduced  from  the  .04cm  mesh  to  the  ,01cm  mesh  but  still 
readily  visible.  The  mesh  was  refined  one  more  time  to  .005cm  to  observe  whether  this 
banding  could  be  reduced  further.  The  pressure  contour  at  50  microseconds  is  shown  in 
Figure  4.19  for  the  .005cm  mesh  to  allow  comparison  with  the  previous  data.  The 
amount  of  material  banding  is  slightly  reduced  again  but  still  readily  apparent.  To  refine 
the  mesh  any  further  would  require  considerable  computational  time  that  is  not  currently 
feasible  since  the  Boundary  Layer  method  is  limited  to  single  processor  computing. 
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Normal  Impact  Results,  Part  2 

The  results  obtained  for  a  normal  impact  at  1  km/s  with  a  Vascomax  300 
projectile  and  a  1080  Steel  target  reveal  an  event  that  is  highly  dominated  by  plastic 
deformation.  Pressures  near  the  impact  surface  exceed  2  GPa,  the  uniaxial  compressive 
yield  stress  of  the  harder  projectile  material.  The  impact  sequence  for  the  No  Slide  Line 
case  is  shown  in  Figures  4.20  and  4.21  by  way  of  pressure  contours  at  10  ps  intervals. 
The  limited  penetration  and  severe  deformation  of  the  projectile  is  attributed  to  the 
similarity  of  the  material  properties  which  result  in  similar  material  responses.  However, 
greater  defonnation  is  observed  in  the  target  material  since  the  yield  strength  is  just  over 
1/3  that  of  the  projectile.  By  80  ps,  as  it  comes  to  rest,  the  amount  of  projectile 
deformation  is  quite  extensive.  There  is  also  an  extensive  zone  of  negative  pressure  that 
has  developed  in  the  projectile.  Just  visible  at  10  ps  is  a  band  of  high  pressure  that 
develops  in  the  projectile  which  increases  in  size  and  magnitude.  This  high  pressure 
band  propagates  towards  the  top  of  the  projectile  as  the  projectile  continues  to  defonn 
and  mushroom  radially  outwards.  The  formation  of  the  negative  pressure  region  begins 
at  approximately  50  ps,  just  below  the  high  pressure  band,  and  expands  as  the  high 
pressure  band  propagates  upwards.  It  is  hypothesized  that  the  negative  pressure  region 
forms  due  to  a  pressure  release  effect  cause  by  the  propagation  of  the  high  pressure  band 
combined  with  an  amplified  elastic  response  due  to  the  similarity  of  the  two  materials 
involved.  The  fact  that  this  negative  pressure  region  is  most  extensive  in  the  No  Slide 
Line  case  supports  the  hypothesis  that  some  material  coupling  effect  is  occurring. 
Whether  this  level  of  plastic  deformation  would  be  observed  in  an  experimental  impact 
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test  is  highly  questionable.  This  extensive  plastic  deformation  is  most  likely  a  purely 
numerical  construct,  since  a  fracture  criteria  was  not  included  in  these  simulations. 

Similar  deformation  is  observed  for  the  Slide  Line  and  Boundary  Layer  cases  as 
seen  in  Figures  4.22  -  4.25.  The  only  significant  difference  observed  between  the  three 
cases  is  the  rate  at  which  the  nose  of  the  projectiles  mushrooms.  The  deformation  of  the 
projectile  for  the  Boundary  Layer  case  appears  to  develop  at  the  slowest  rate.  This  is 
most  visible  at  30  ps  where  a  gap  exists  between  the  advancing  front  of  the  projectile 
material  and  the  trailing  face  of  the  deformed  target  material.  However,  this  gap  does  not 
persist  as  the  projectile  continues  to  mushroom  and  reestablish  contact  with  the  target 
material.  In  all  three  cases,  the  projectile  is  observed  to  contact  the  target  as  the  impact 
event  comes  to  an  end.  The  difference  in  deformation  rate  is  also  evident  at  60  ps  by 
comparing  the  deformed  height  of  the  projectile.  While  the  penetration  depth  is 
essentially  identical,  the  undeformed  aft  end  of  the  projectile  is  right  at  the  2cm  mark  in 
the  Boundary  Layer  case,  while  it  is  clearly  below  the  2cm  mark  in  the  other  two  cases. 
These  two  observations  suggest  that  the  boundary  layer  algorithm  is  artificially 
increasing  the  yield  stress  of  the  projectile.  Since  it  was  implemented  to  preserve  the 
penetrator  shape  in  low  velocity  impacts,  perhaps  that  effect  is  appearing  in  this  case 
inappropriately.  Perhaps  this  technique  is  not  appropriate  for  use  when  the  projectile  is 
known  to  undergo  substantial  deformation. 
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Figure  4.20  Pressure  Contours  at  10,  20,  30,  and  40  jjs,  No  Slide  Line 
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Figure  4.21  Pressure  Contours  at  50,  60,  70,  and  80  gs,  No  Slide  Line 
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Figure  4.23  Pressure  Contours  at  50,  60,  70,  and  80  ps,  Slide  Line 


51 


Pressure 


Pressure 


2DC 

Vasco300  ->  1080, 


X  (cm) 

1000  m/s,  mu  =  0.2,  dx=0.01 


2/11/05  18:19:18  CTH 


Tim  6= 1.00  04x10“ 


2DC 

Vasco3Q0  ->  1080,  1000  m/s, 
BKRCYU  2/11/05  19:30:20 


X  (cm) 

mu  =  0.2,  dx=0.01 
CTH  2580 


n  6= 2. 0005  7x10“ 


2DC  X  (cm) 

Vasco300  ->  1080,  1000  m/s,  mu  =  0.2,  dx=0.01 
BKRCYU  2/11/05  20:41:25  CTH  3839 


Tim6=3. 0008x10” 


2DC  X  (cm) 

Vasco3Q0  ->  1080,  1000  m/s,  mu=0.2,  dx=0.01 

BKRCYU  2/11/05  21:53:23  CTH  5089  TJm6=4.00006x10“ 


Figure  4.24  Pressure  Contours  at  10,  20,  30,  and  40  jjs,  Boundary  Layer 
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Figure  4.25  Pressure  Contours  at  50,  60,  70,  and  80  ps,  Boundary  Layer 
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The  assertion  that  the  extensive  defonnation  of  both  projectile  and  target  material 
is  primarily  plastic  deformation  is  supported  by  the  temperature  profiles  shown  in  Figures 
4.26  -  4.28.  The  temperature  rise  in  the  majority  of  the  defonned  material  never  reaches 
the  melting  temperatures  of  1650  K  for  the  projectile  and  1835  K  for  the  target.  The  No 
Slide  Line  case  does  show  a  high  temperature  zone  along  the  entire  contact  surface  of  the 
projectile  which  is  not  present  in  the  other  cases.  This  is  consistent  with  the  restriction 
imposed  by  this  technique  that  the  materials  are  joined  at  the  contact  surface.  The  same 
level  of  plastic  strain  in  the  cells  at  the  contact  surface  would  result  in  a  higher 
temperature  in  the  projectile  material  due  to  its  higher  yield  strength.  Consideration  of 
friction  mechanisms  identified  by  Bayer[3]  would  suggest  that  the  No  Slide  Line 
technique  simulates  the  extreme  scenario  of  pure  adhesion  along  perfectly  smooth 
surfaces  with  100%  contact.  Likewise,  the  Boundary  Layer  technique  is  simulating 
adhesion  along  surfaces  with  less  than  100%  contact;  while  the  Slide  Line  technique  as 
intended  is  simulating  frictionless  sliding  to  the  extent  possible  considering  numerical 
artifacts  may  be  a  factor. 
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Figure  4.26  Temperature  at  30  ps,  No  Slide  Line 
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Figure  4.27  Temperature  at  30  ps,  Slide  Line 
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Figure  4.28  Temperature  at  30  ps,  Boundary  Layer 


To  confirm  the  similarities  between  the  three  material  interface  methods  that  have 
been  observed  on  a  global  scale,  two  projectile  material  tracers  were  selected  to  evaluate 
the  conditions  in  the  projectile  near  the  contact  surface.  These  two  tracers  are  numbered 
86,  with  initial  coordinates  x=.205  and  y=.305,  and  107  with  initial  coordinates  x=.355 
and  y=.505.  Tracer  86  is  representative  of  the  area  near  the  center  of  the  projectile  nose 
while  tracer  107  is  representative  of  the  area  near  the  outer  edge  of  the  nose  which  sees 
the  greatest  deformation.  These  tracers  are  located  near  enough  to  the  contact  surface  to 
evaluate  the  influence  of  the  different  material  interface  methods  on  the  solution  but  far 
enough  to  not  be  directly  involved  in  the  unique  numerical  processes  employed  by  each 
method. 
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Before  discussing  these  figures,  the  issue  of  thermodynamic  warnings  mentioned 
previously  needs  to  be  addressed  again.  The  solution  obtained  with  the  Slide  Line 
method  automatically  aborted  at  approximately  95  ps  due  to  the  excessive  number  of 
thennodynamic  warnings,  exceeding  the  arbitrary  limit  of  one  million.  There  is  a  note  in 
the  CTH  reference  manual  that  suggests  caution  when  using  the  Slide  Line  method  since 
it  can  generate  unexpected  results.  With  this  in  mind  the  results  presented  for  the  Slide 
Line  method  should  be  viewed  with  some  concern  for  their  accuracy. 

The  displacement  history  of  tracer  86  (Figure  4.29)  shows  a  similar  response 
pattern  for  all  three  methods  but  a  significant  difference  in  magnitude.  The  largest  final 
displacement  is  observed  with  the  No  Slide  Line  method  (.37cm),  followed  by  the  Slide 
Line  method  (.32cm),  and  the  least  displacement  is  observed  with  the  Boundary  Layer 
method  (.22cm).  This  supports  the  earlier  observation  that  less  plastic  deformation 
occurs  with  the  Boundary  Layer  method.  The  largest  difference  in  final  displacement  is 
about  .15cm  or  40%  of  the  largest  displacement  observed.  A  small  recovery  is  also 
evident  once  the  maximum  displacement  has  been  reach  with  the  Slide  Line  and 
Boundary  Layer  methods. 

Similar  observations  can  be  made  regarding  the  temperature  history  of  tracer  86 
(Figure  4.30).  The  percent  difference  in  temperature  increase  is  not  as  large  however 
across  the  three  methods.  The  difference  of  about  125K  is  only  15%  of  the  maximum 
temperature  increase. 

The  pressure  history  of  tracer  86  (Figure  4.3 1)  indicates  an  initial  spike  in 
pressure  due  to  the  initial  impact  which  relaxes  until  another  pressure  rise  is  observed  at 
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40  ps.  While  the  initial  spike  has  a  similar  magnitude  of  about  2GPa  across  all  three 
methods,  the  magnitude  resulting  from  the  second  pressure  rise  varies.  The  peak 
magnitude  of  this  second  pressure  rise  is  approximately  2GPa  for  the  No  Slide  Line 
method,  3GPa  for  the  Slide  Line  method  and  3.3GPa  for  the  Boundary  Layer  method. 
This  second  pressure  rise  is  possibly  due  to  a  pressure  wave  generated  due  to  the  initial 
impact  which  has  propagated  thru  the  projectile  and  reflected  back  from  the  free  surface 
at  the  tail  end.  The  significant  differences  in  the  magnitude  of  this  reflected  pressure 
wave  is  unexpected. 

The  deviatoric  shear  stress  history  of  tracer  86  (Figure  4.32)  also  shows  similar 
responses  for  all  three  methods.  An  initial  negative  shear  spike  of  about  -900MPa  due  to 
the  impact  is  observed  followed  by  a  rapid  increase  in  shear  to  an  equally  large  positive 
value.  The  positive  shear  is  sustained  until  approximately  30  ps,  when  the  shear  values 
rapidly  drop  and  remain  negative  for  the  remainder  of  the  solution  time.  This  drop  results 
in  shear  magnitudes  less  than  half  that  of  the  initial  spike  for  the  Boundary  Layer  method 
while  the  drop  is  more  severe  for  the  other  two  methods  and  is  approximately  the  same 
magnitude  as  the  initial  spike.  This  pattern  of  shear  stress  development  could  be  a  result 
of  the  differences  in  the  yield  stress  and  plastic  deformation  rates  of  the  two  materials.  A 
shear  response  would  be  expected  if  the  materials  which  are  in  contact  are  deforming  at 
different  rates. 

The  y-component  of  nonnal  stress  for  tracer  86  (Figure  4.33)  shows  significant 
time  response  differences  for  the  three  methods  from  10  to  60  ps.  With  the  No  Slide  Line 
method  the  response  is  almost  identical  to  the  pressure  response  shown  in  Figure  4.3 1 .  A 
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slight  rise  is  evident  with  the  Slide  Line  method  that  follows  the  pressure  response  but 
not  to  the  extent  of  the  No  Slide  Line  method.  The  Boundary  Layer  method  shows 
relatively  little  change  during  this  time  interval.  After  60  ps,  the  Slide  Line  method  and 
the  Boundary  Layer  method  show  a  steady  drop  toward  zero  stress  while  the  No  Slide 
Line  method  only  drops  to  about  lGPa. 
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Figure  4.29  Tracer  86,  x-Coordinate  History 
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TEMPERATURE  (1Q5K) 
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Figure  4.30  Tracer  86,  Temperature  History 
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PRESSURE  (10SD/CM2) 
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Figure  4.31  Tracer  86,  Pressure  History 
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YYSTRESS  (10J[)/CM2) 
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The  displacement  plots  of  tracer  107  (Figure  4.34)  show  displacement  magnitudes 
that  are  very  similar  for  all  three  methods  with  final  displacements  differing  by  only  .02cm  or 
about  6%  of  the  greatest  displacement  observed.  The  general  displacement  response  is 
similar  to  that  observed  with  tracer  86  except  for  an  interruption  observed  in  the  20  to  40  ps 
time  interval  for  the  No  Slide  Line  method.  Also,  the  small  recovery  after  maximum 
displacement  was  reached  that  was  observed  previously  is  almost  nonexistent  for  this  tracer. 

The  temperature  history  for  tracer  107  (Figure  4.35)  also  shows  overall  temperature 
responses  similar  to  those  observed  for  tracer  86  with  the  exception  again  of  the  No  Slide 
Line  method.  A  disruption  in  the  temperature  increase  rate  in  the  same  20  to  40  ps  time 
interval  for  the  No  Slide  Line  method  is  observed.  The  temperature  becomes  very  erratic  in 
this  time  interval  but  still  results  in  a  steady  final  temperature  of  approximately  1225K.  A 
much  smaller  disruption  in  the  temperature  development  for  the  Slide  Line  method  is  also 
observed  over  the  same  time  interval  and  the  final  temperature  is  approximately  975K,  the 
smallest  temperature  increase.  The  temperature  development  for  the  Boundary  Layer  method 
appears  relatively  smooth  and  reaches  a  final  temperature  of  approximately  1 100K.  While  a 
temperature  difference  of  only  125K  was  observed  between  the  three  methods  for  tracer  86,  a 
higher  difference  of  250K,  which  is  27%  of  the  largest  temperature  increase,  is  observed  for 
tracer  107. 

Figure  4.36  shows  three  very  distinct  pressure  histories  for  tracer  107.  Negative 
pressure  values  are  indicated  in  the  20-40  ps  time  interval  for  the  No  Slide  Line  method 
while  pressure  values  remain  clearly  positive  for  the  other  two  methods.  The  maximum 
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pressure  of  2.5GPa  occurs  with  the  Slide  Line  method  but  all  three  methods  do  show  a  trend 
towards  a  final  pressure  value  of  lGPa. 

The  deviatoric  shear  stress  history  of  tracer  107  (Figure  4.37)  is  roughly  similar  to 
that  of  tracer  86  (Figure  4.32)  although  no  initial  negative  shear  spike  is  observed. 
Additionally,  compared  to  the  shear  stress  magnitudes  observed  for  tracer  86,  the  final  shear 
stress  magnitude  is  twice  as  large  for  the  Boundary  Layer  method,  but  less  than  one  half  of 
the  magnitude  for  the  Slide  Line  method,  and  only  one  fourth  of  the  magnitude  for  the  No 
Slide  Line  method. 

The  y-component  of  nonnal  stress  for  tracer  107  (Figure  4.38)  shows  a  time  response 
very  similar  to  the  pressure  response.  There  is  also  more  similarity  between  results  for  the 
Slide  Line  method  and  the  Boundary  Layer  method.  The  largest  peak  value  of  3.5GPa 
occurs  with  the  Slide  Line  method  but  the  No  Slide  Line  method  results  in  the  largest  final 
value  of  1.8GPa.  It  should  also  be  noted  that  these  stress  values  are  consistently  greater  in 
magnitude  than  the  pressure  values. 

These  tracer  particle  history  plots  would  suggest  that  while  the  three  different 
material  interface  methods  do  not  produce  significantly  different  overall  effects  in  an  impact 
scenario,  there  are  significant  differences  that  exist  on  a  smaller,  localized  scale. 

The  reason  for  the  curious  results  observed  with  the  No  Slide  Line  method  in  the  20  - 
40  ps  time  interval  is  not  readily  apparent  but  looking  back  at  Figure  4.26  reveals  a  projectile 
deformation  state  that  may  suggest  a  cause.  Since  the  two  materials  are  joined  at  the  contact 
surface  in  the  No  Slide  Line  method,  projectile  material  is  pulled  along  with  the  target 
material  and  folds  over  on  itself  in  the  20-40  ps  time  interval.  This  secondary  impact  would 
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definitely  affect  the  displacement  rate  of  material  in  the  vicinity  of  tracer  107.  Since  the  two 
other  methods  allow  for  relative  material  sliding,  this  effect  is  not  observed. 
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TEMPERATURE  (105K) 
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PRESSURE  (10SD/CM2) 
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Pure  Sliding  Results 

The  scenario  of  pure  sliding  reveals  severe  difficulties  with  the  Slide  Line 
method.  Since  this  method  establishes  a  frictionless  sliding  interface  at  the  contact 
surface,  we  would  not  expect  any  shear  stress  to  develop.  However,  a  plot  of  the 
deviatoric  shear  stress  reveals  the  development  of  a  shear  response  throughout  the 
projectile  and  target  within  the  first  few  microseconds  as  can  be  seen  in  Figure  4.39.  This 
development  must  be  attributed  to  numerical  noise.  As  the  solution  develops,  this 
numerical  noise  dominates  the  solution  and  begins  to  generate  completely  nonphysical 
results.  Thermodynamic  warnings  are  prevalent  again  and  the  solution  actually 
tenninates  at  approximately  49  ps  when  the  number  of  warnings  exceeds  one  million. 

The  dark  blue  cells  seen  in  Figure  4.40  of  the  temperature  profile  indicates  regions  where 
the  temperature  has  dropped  below  ambient.  With  no  heat  transfer  considered,  this  must 
surely  be  a  nonphysical  result.  Aside  from  the  fact  that  the  Slide  Line  method  is 
inappropriate  for  this  particular  scenario,  any  other  comments  would  be  unwarranted. 
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Plots  of  the  deviatoric  shear  stress  for  the  other  two  methods  reveal  two  very 
distinct  developments.  Figures  4.41a  and  4.41b  show  shear  stress  states  at  24 
microseconds  and  50  microseconds  for  the  No  Slide  Line  method.  The  shear  response 
seen  as  alternating  vertical  regions  of  positive  and  negative  shear  is  extensive  in  both  the 
projectile  and  the  target.  As  the  solution  progresses,  this  shear  pattern  in  the  target  and 
the  shear  response  itself  weakens  and  begins  to  disappear,  while  the  shear  pattern  in  the 
projectile  appears  to  intensify.  In  addition,  the  shear  response  in  the  projectile  also  shows 
alternating  horizontal  regions  of  positive  and  negative  shear  originating  near  the 
projectile  nose.  The  shear  development  for  the  Boundary  Layer  method  shown  in  figures 
4.42a  and  4.42b  show  a  very  intense  region  of  shear  in  the  target  material  near  the 
material  interface.  This  should  be  the  effect  of  the  friction  body  force  which  is  added 
through  the  Boundary  Layer  Algorithm.  The  lack  of  other  significant  shear  response 
features  is  attributed  to  the  slip  layer  which  would  not  transfer  shear  into  the  target  and 
produce  the  material  response  seen  with  the  No  Slide  Line  method. 

The  shear  response  pattern  could  also  be  a  result  of  numerical  noise  due  to  mesh 
size  issues.  The  mesh  cell  size  of  .01cm  for  this  case  was  carried  over  from  the  impact 
studies  and  may  not  be  sufficient  to  capture  the  true  material  response  for  this  scenario. 
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Figure  4.41a  Deviatoric  Shear  Stress  at  24  gs,  No  Slide  Line 
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dx=0.01 

Time= 5.00052x10" 


Figure  4.41b  Deviatoric  Shear  Stress  at  50  gs,  No  Slide  Line 
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Stress  Deviator 
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mu=0.2, 

3084 


0.8  1  . 

dx  =  0.01 

Tim  e=2. 40021x10" 


Figure  4.42a  Deviatoric  Shear  Stress  at  24  gs,  BOUNDARY  LAYER 
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1  .  □ 


X  (cm) 

Vasco300  Sliding  Inside  1080,  1000  m/s,  mu=0.2,  dx  =  0.Q1 

BQRAKA  2/17/05  19:11:22  CTH  6343  Time=5.00054x10" 5 


Figure  4.42b  Deviatoric  Shear  Stress  at  50  gs,  BOUNDARY  LAYER 


77 


The  temperature  development  for  the  No  Slide  Line  method  and  the  Boundary 
Layer  method  are  similar  in  magnitude  and  distribution  as  shown  in  Figures  4.43a  and 
4.43b.  While  the  No  Slide  Line  method  produces  a  smoother  temperature  distribution, 
both  methods  produce  a  high  temperature  projectile  surface  layer  near  the  leading  edge 
which  transitions  into  the  target  surface.  It  is  somewhat  surprising  that  the  two 
fundamentally  different  methods  would  produce  such  similar  results.  The  dissipation  of 
shear  stress  in  the  target  observed  with  the  No  Slide  Line  method  along  with  the 
established  layer  of  high  temperature  target  material  suggests  the  formation  of  a  fluid 
layer.  As  this  fluid  layer  grows,  less  shear  stress  is  transferred  into  the  target.  It  would 
appear  that  given  enough  time,  both  the  No  Slide  Line  method  and  Boundary  Layer 
method  will  generate  the  same  solution. 


T  emperature 


0  .  0 


2DC  X  (cm) 

Vasco3Q0  Sliding  Inside  1080,  1000  m/s,  noSL, 
BQTBYR  2/17/05  21:19:45  CTH  6331 


Temperature 

GO 
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1.613x10* 
1.S26xl03 
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1.351xl03 
1.264-xl  03 
1.177x10s 
1.090x10s 
1.003x10s 
9.154x10s 
8.282x10s 
7.410x10s 
6.538x10  s 
5.667x10s 
_  4.795x10s 

1  3.923x10s 
3.051x10s 
2.179x10s 
1.306x10s 
4.359x10' 


dx=0.01 

Time=5.00052x10-5 


Figure  4.43a  Temperature  at  50  ps,  No  Slide  Line 
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T  emperature 
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Vasco300  Sliding  Inside  1080,  1000 
BQRAKA  2/17/05  19:11:22  CTH 


X  (cm) 
m/s,  mu=0.2, 
6343 


dx  =  0.01 

Time=5.00054x10"  5 


Figure  4.43b  Temperature  at  50  ps,  Boundary  Layer 
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V.  Conclusions  and  Recommendations 


Chapter  Overview 

Through  the  course  of  this  thesis  research,  computational  methods  for  the 
treatment  of  material  interfaces  available  in  the  hydrocode  CTH  were  investigated  with 
the  goal  of  establishing  some  criteria  to  aid  in  the  application  of  these  methods  to 
scenarios  for  which  they  were  not  originally  intended.  The  three  methods  investigated 
are  repeated  here: 

1.  No  Slide  Line,  materials  are  joined  at  the  contact  surface 

2.  Slide  Line,  a  frictionless  sliding  interface  is  implemented  at  the  contact  surface 

3.  Boundary  Layer,  friction  is  included  and  the  sliding  interface  is  shifted  into  the 
softer  material 

The  specific  application  of  interest  is  the  simulation  of  the  hypervelocity  gouging 
phenomenon  observed  at  the  HHSTT. 

Conclusions  of  Research 

The  first  phase  of  this  research  analyzed  a  deep  penetration  impact  event 
involving  a  steel  projectile  and  an  aluminum  target.  The  projectile  had  a  yield  stress 
approximately  six  times  that  of  the  target  and  had  nearly  3  times  the  stiffness.  The 
resulting  impact  clearly  showed  an  event  dominated  by  target  plastic  deformation.  For  a 
relatively  coarse  mesh  of  .04cm,  these  impact  studies  indicated  that  the  Boundary  Layer 
method  produced  results  that  were  the  closest  to  experimental  data.  However,  after 
refining  the  mesh  to  .01cm,  the  results  were  almost  identical  for  all  three  methods.  The 
effect  of  material  interface  phenomena  in  this  scenario  did  not  appear  to  be  a  significant 
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contributor  to  the  impact  solution.  For  this  type  of  scenario  then,  the  suggestion  would 
be  to  use  the  No  Slide  Line  method  since  it  takes  no  effort  on  the  part  of  the  user  to 
implement  and  requires  the  least  computational  effort. 

The  second  phase  of  this  research  analyzed  a  normal  impact  event  involving  a 
Vascomax300  projectile  and  a  1080  steel  target,  in  an  attempt  to  model  materials  in  use  at 
the  HHSTT.  The  projectile  had  a  yield  stress  approximately  three  times  that  of  the  target 
but  a  stiffness  value  only  15%  greater.  The  resulting  impact  showed  significant  plastic 
deformation  of  both  the  projectile  and  the  target.  Material  interface  effects  were 
observed  to  produce  detectable  differences  in  the  CTH  solution.  The  No  Slide  Line 
method  produced  excessive  defonnation  of  the  projectile  in  this  scenario.  Since  the 
projectile  was  artificially  joined  to  the  target  material  at  their  contact  surfaces,  the  higher 
plastic  defonnation  rate  of  the  target  material  pulled  projectile  material  along  with  it. 

The  Slide  Line  method  produced  results  similar  to  the  Boundary  Layer  method  although 
the  issue  of  thermodynamic  warnings  forced  this  solution  to  abort  early  which  raises 
concerns  about  the  accuracy  of  that  solution.  The  Boundary  Layer  method  produced 
results  that  were  consistent  with  the  other  methods  from  an  overall  perspective  but  did 
not  show  any  significant  erratic  phenomena  in  the  tracer  histories  that  were  present  in  the 
other  two  methods.  Even  though  it  was  observed  that  the  Boundary  Layer  method  was 
artificially  increasing  the  yield  strength  of  the  projectile,  the  effect  was  relatively  small. 
The  Boundary  Layer  method  appears  to  be  the  method  of  choice  in  this  scenario  since  it 
is  numerically  stable  and  the  user  has  some  options  with  this  method  to  fine  tune  the 
solution  to  reflect  experimental  data  if  available.  Once  a  sufficiently  similar  simulation  is 
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established,  other  aspects  of  the  event  of  interest  which  were  not  physically  observable 
can  be  studied  in  simulations  with  greater  confidence. 

The  third  phase  of  this  research  analyzed  a  sliding  scenario  involving  a 
Vascomax300  projectile  and  1080  steel  target.  This  extreme  scenario  which  attempts  to 
isolate  the  ability  of  the  three  methods  to  model  pure  sliding  produced  some  unexpected 
results.  The  numerical  instability  of  the  Slide  Line  method  was  made  quite  evident  with 
this  scenario.  The  similarities  between  the  No  Slide  Line  method  and  the  Boundary 
Layer  method  were  not  at  ah  expected.  The  issue  of  mesh  refinement  however  precludes 
any  definitive  comparison  of  the  two  methods.  Although  the  results  do  suggest  that  if  the 
relative  tangential  velocity  between  two  materials  is  sufficiently  high,  then  both  methods 
are  suitable  for  use  which  would  lead  to  the  recommendation  of  the  No  Slide  Line 
method  simply  for  ease  of  implementation  and  reduction  in  computational  effort. 

Based  on  these  limited  findings,  the  No  Slide  Line  method  would  be 
recommended  for  general  use  in  high  energy  impact  studies.  Due  to  the  numerical 
instability  encountered  with  the  Slide  Line  method,  its  application  should  be  limited.  The 
Boundary  Layer  method  is  too  limited  in  its  application  for  2D  models.  The  lack  of 
support  for  parallel  data  processing  is  also  a  major  drawback  with  this  method. 
Implementing  it  in  a  2D  non-axisymmetric  problem  won’t  generate  an  error  and  cause  the 
CTH  run  to  abort  but  the  solution  will  not  include  any  friction  input  the  user  may  be 
interested  in.  Unless  the  problem  is  2D  axisymmetric  and  results  with  a  relatively  coarse 
mesh  are  acceptable,  so  that  solutions  can  be  obtained  with  a  standalone  workstation,  the 
Boundary  Layer  method  offers  no  advantage. 
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This  leads  to  a  recommendation  for  simulations  of  the  hypervelocity  gouging 
problem  at  the  HHSTT.  For  the  2D  plane  strain  model  currently  in  use,  the  No  Slide 
Line  method  should  be  implemented.  If  there  is  interest  in  applying  the  Boundary  Layer 
method  for  its  capability  to  explicitly  add  a  friction  force  to  the  problem  then  a  3D  model 
will  need  to  be  created. 

Recommendations  for  Future  Research 

This  research  has  by  no  means  been  exhaustive  in  its  attempt  to  discern  which 
material  interface  methods  are  suitable  for  which  situations.  The  recommendations  for 
future  research  are: 

1 .  Evaluate  the  three  methods  in  oblique  impact  scenarios. 

2.  Evaluate  the  three  methods  in  impact  scenarios  with  coatings. 

3.  Evaluate  the  three  methods  for  a  3D  impact  model. 

Any  future  research  effort  should  also  implement  the  heat  transfer  capabilities  of  CTH  in 
order  to  capture  the  extent  to  which  friction  affects  the  thermodynamic  state  of  the 
materials. 
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Appendix:  CTH  Input  File  for  Impact  Scenario  Two 


*eor*  cgenin 
* 

*This  is  the  combined  CTHGen/CTH  input  file 
*This  is  the  CTHGen  input  section 
* 

Vasco300  ->  1080,  1000  m/s,  mu=0.2,  dx=0.01 
* 

*A11  output  from  CTHGen  will  have  this  title 

*This  is  the  impact  scenario  of  a  Vascomax300  projectile  impacting  a 
*1080  steel  target  at  1000  m/s,  implementing  the  boundary  layer  method 
*with  a  coefficient  of  friction  of  0.2,  and  a  .01cm  mesh  cell  size  in 
*the  area  of  interest 
* 

control 

mmp 

ep 

vpsave 

endcontrol 

* 

mesh 

block  geom=2dc  type=e 

* 

*  2D-cylindrical  coordinates,  eulerian  mesh 

*  axisymmetric  problem,  so  only  half  model  generated 

* 

x0=0 . 

xl  dxf=0 . 01  dxl=0 . 01  w= . 5 
x2  dxf=0 . 01  dxl=0 . 04  w=5 . 5 
endx 

* 

*x  ranges  from  0.  to  6. 

*dxf  is  size  of  first  cell  in  the  range 

*dxl  is  size  of  last  cell  in  the  range 
* 

y0=-10 . 0 

yl  dyf=0.04  dyl=0.04  w=4 . 
y2  dyf=0 . 04  dyl=0.01  w=2 . 
y3  dyf=0 . 01  dyl=0.01  w=6. 
y4  dyf=0 . 01  dyl=0.04  w=2 . 
y4  dyf=0 . 04  dyl=0.04  w=3 . 5 
endy 

* 

*y  ranges  from  -10.  to  7.5 

*dyf  is  size  of  first  cell  in  the  range 

*dyl  is  size  of  last  cell  in  the  range 
* 

xactive  0  0.5 
yactive  -1  7.5 

* 

*this  is  the  initial  mesh  region  containing  activity 

*helps  speed  up  the  computation  by  not  working  on  the  entire  mesh  from 
*the  start 
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endb 

endmesh 


insertion  of  material 
package  projectile 

* 

*  main  body  of  Vasacomax300  projectile,  7.11mm  diameter,  71.1mm  length 

* 

material  2 
numsub  50 
xvel  0.0 
yvel  -1000. e2 

* 

*  cgs  units,  cm/s 

* 

insert  box 

pi  0.  0.3555  p2  0.3555  7.1476 
endinsert 

* 

*  only  one  insert  allowed  per  package 

* 

endpackage 
package  projectile 

* 

*  hemispherical  nose  of  the  Vasacomax300  projectile 

* 

material  2 
numsub  100 
xvel  0.0 
yvel  -1000. e2 
insert  circle 
cen  0.  0.3555  rad  0.3555 
endinsert 
endpackage 
package  target 

* 

*  1080  Steel  target,  152mm  diameter,  198mm  length 

* 

material  1 
numsub  100 
xvel  0. 
yvel  0. 
insert  box 
pi  0.  -19.8 
p2  7.6  0 
endinsert 
endpackage 
endinsertion 
* 

edit 

expanded 

endedit 
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eos 

matl  ses  iron 
mat2  ses  Steel_V300 

* 

*  Sesame  Equation  of  State  tables,  limited  material  selection 

* 

endeos 

* 

epdata 
vpsave 
mix  3 

* 

*calculate  yield  strength  in  mixed  cells  using  volume 
*averaged  yield  strength  normalized  by  the  sum  of  the 
*volume  fractions  of  the  materials  that  can  support  shear 
* 

matep  1 

johnson-cook  IRON 
poisson  0.28 
yield=7 . 0e9 

* 

*Viscoplastic  rate  dependent  model  of  Iron,  closest  material  available 
*to  1080  steel 
* 

matep  2 
st=l  9 

poisson=0 . 27 
yield=14 . 47e9 

* 

*Viscoplastic  rate  independent  model  of  Vascomax 
*only  some  Steinberg-Guinan-Lund  models  are  rate  dependent 
* 

BLINT  1  soft  1  hard  2  csl  1.5  cbl  1.5  fric  0.2  corr 

* 

*Boundary  Layer  Interface  model 

*first  number  is  an  identifier  for  each  BLINT  entry,  if  you  have  3 

*materials,  you  will  need  two  BLINT  entries  uniquely  identified 

*soft  1,  material  identifier  of  the  soft  material 

*hard  2,  material  identifier  of  the  hard  material 

*csl  1.5,  slip  layer  width  is  1.5  cell  diagonals 

*cbl  1.5,  boundary  layer  width  is  1.5  cell  diagonals 

*fric  0.2,  coefficient  of  friction  is  0.2 

*corr,  increases  hard  material  yield  stress  by  a  factor  proportional  to 

*  the  increased  radius  of  the  projectile  due  to  soft  material  gain 

*  in  mixed  material  cells 

* 

ende 

* 

tracer 

* 

*tracer  input,  starting  x,y  coord,  to  ending  x,y  coord. 

*n  =  total  number  of  tracers  to  distribute  including  the  end  pts. 

* 

*Centerline 
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add  0.  0.  to  0.  7.0  n  8 
add  0.  7.25 


*projectile  nose  boundary  layer 


add 

.005 

.005 

to 

.405 

.005 

n=9 

add 

.005 

.015 

to 

.405 

.015 

n=9 

add 

.005 

.025 

to 

.405 

.025 

n=9 

add 

.005 

.035 

to 

.405 

.035 

n=9 

add 

.005 

.045 

to 

.405 

.045 

n=9 

add 

.005 

.055 

to 

.405 

.045 

n=9 

add 

.005 

.  105 

to 

.405 

.  105 

n=9 

add 

.005 

.205 

to 

.405 

.205 

n=9 

add 

.005 

.305 

to 

.405 

.305 

n=9 

add 

.005 

.405 

to 

.405 

.405 

n=9 

add 

.005 

.505 

to 

.405 

.505 

n=9 

*target  boundary  layer 


add 

.005 

-.005 

to 

.405 

-.005 

n=9 

add 

.005 

-.015 

to 

.405 

-.015 

n=9 

add 

.005 

-.025 

to 

.405 

-.025 

n=9 

add 

.005 

-.035 

to 

.405 

-.035 

n=9 

add 

.005 

-.045 

to 

.405 

-.045 

n=9 

add 

.005 

-.055 

to 

.405 

-.055 

n=9 

add 

.005 

-.105 

to 

.405 

-.105 

n=9 

add 

.005 

-.205 

to 

.405 

-.205 

n=9 

add 

.005 

-.305 

to 

.405 

-.305 

n=9 

add 

.005 

-.405 

to 

.405 

-.405 

n=9 

add 

endt 

.005 

-.505 

to 

.405 

-.505 

n=9 

matplot 
x=0  . ,  .  5 

y=-l . , 1 . 

Endm 

* 

*generates  a  plot  of  materials  in  the  metamat  file 

*use  popxll  command  to  view  file  or  popjpg  command  to  make  a  jpeg 

* 

matplot 
x=0 . , 5 . 
y=~3 . ,  7 
endm 

'k'k'k'k'k'k'k-k'k'k'k'k'k'k'k-k'k'k'k-k'k'k'k-k'k'k'k'k'k'k'k-k'k'k'k-k'k-k'k-k'k'k'k'k'k'k'k-k'k'k'k-k'k'k'k-k'k-k'k'k 

*eor*  cthin 
* 

*CTH  input  section 
* 

Vasco300  ->  1080,  1000  m/s,  mu=0.2,  dx=0.01 
* 

*A11  output  from  CTH  will  have  this  title 
* 

control 

mmp 

tstop  =  200.0e-6 
nscycle  =  150000 
cpshif t=600 . 
ntbad=1000000 
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*mmp  allows  multiple  pressures  and  temperatures  in  mixed  material  cells 
*tstop  =  solution  end  time 

*nscycle  =  solution  will  end  if  number  of  cycles  is  exceeded 
*cpshift  =  allows  extra  time  for  CTH  to  write  data 
*ntbad  =  maximum  thermodynamic  warnings  before  CTH  aborts 
* 

endc 

* 

convct 

interf ace=high  resolution 
endc 
* 

edit 

shortt 

time=0.  dt=20.e-6 

* 

*restart  data  will  be  written  every  "dt"  seconds 
* 

ends 

longt 

tim=0  dt=l . 0 

* 

*restart  data  will  be  written  every  "dt"  seconds 
* 

endl 

histt 

time=0  dt=0.5e-6 
htracer  all 

* 

*history  data  will  be  written  to  hcth  every  "dt"  seconds 
*data  for  all  tracer  points  is  written 
* 

endh 

plott 

time=0  dt=l.e-6 
time=50.e-6  dt=2.e-6 
time=100.e-6  dt=5.e-6 

* 

*cthplot  data  is  written  to  restart  file  every  "dt"  seconds 
*starting  at  "time" 

*beware  restart  file  size  limit  is  2GB 
endp 
ende 
* 

boundary 

bhydro 

bxb=0.  bxt=l .  byb=l .  byt=l . 

* 

*symmetry  boundary  condition  at  minimum  x 

*sound  speed  absorbing  boundary  condition  at  max  x,  min  y,  and  max  y 

*equivalent  to  semi-infinite  boundary 
* 

endh 
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endb 

* 

mindt 

time  50.0e-6 
dtmin  1 . Oe-9 
endm 
* 

*CTH  will  abort  if  the  time  step  drops  below  dtmin  after 
* 


"time" 
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